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Abstract 


The  dynamic  instability  characteristics  of  a 
laminated  composite  panel  subjected  to  a  transverse  load  is 
studied  in  this  research.  Up  to  cubic  variations  in  the 
thickness  coordinate  are  included  in  the  inplane 
displacement  field,  and  only  the  constant  component  is  kept 
in  the  transverse  displacement.  The  transverse  shear 
strains  retain  only  linear  displacement  terms  and  vary 
parabolical ly  through  the  thickness  vanishing  at  the  top  and 
bottom  surfaces.  The  complete  quadratic  displacement  func¬ 
tions  are  included  in  the  inplane  strains  to  characterize 
the  large  displacements/rotations  response  during  a  snapping 
process.  A  36  degree  of  freedom  shell  element  is  used  to 
obtain  numerical  results. 

The  static  snap  through  load  versus  dis¬ 
placement  curve,  as  well  as  the  critical  collapse  load,  is 
first  examined  by  invoking  the  Riks  technique  along  with  the 
Newton-Raphson  iteration  scheme  at  each  load  increment 
level.  The  beta-m  time  marching  integration  method  is  then 
employed  to  evaluate  a  dynamic  response.  Two  step  loads, 
with  the  step  magnitude  slightly  below  and  above  the 
critical  collapse  load,  are  introduced  in  the  dynamic 
analysis.  Moreover,  the  damping  effect  is  incorporated  to 


yield  a  steady-state  solution.  The  response  resulting  from 
the  dynamic  steady-state  analysis  subjected  to  a  step  load 
matches  the  displacement  on  the  snap  through  load  versus 
displacement  curve. 


FINITE  ELEMENT  INVESTIGATION  INTO  THE  DYNAMIC  INSTABILITY 


CHARACTERISTICS  OF  LAMINATED  COMPOSITE  PANELS 


I.  INTRODUCTION 


Laminated  composite  materials  are  seeing  widespread  use 
in  many  diverse  industries.  One  of  these  is  the  aerospace 
industry  where,  due  to  the  need  to  minimize  weight,  complex 
shell  configurations  are  common  structural  elements. 
Structural  components  made  from  composite  materials 
typically  have  higher  strength  and  stiffness  to  weight 
ratios  than  those  made  from  isotropic  materials.  In 
addition,  composite  materials’  properties  can  be  tailored  to 
meet  specific  design  goals.  Increased  stiffness  and 
strength  are  designed  only  where  needed.  However,  optimized 
structural  systems  are  often  more  susceptible  to 
instabilities  such  as  buckling,  collapse,  or  vibration, 
especially  if  the  composite  is  a  thin  shell  and  the  load  is 
applied  in  the  transverse  direction.  The  purpose  of  this 
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research  is  to  investigate  numerically  the  static  and 
dynamic  response  of  a  laminated  graphite/ epoxy  composite 
cylindrical  shell  (curved  panel)  subjected  to  transverse 
loading. 

Analytical  Method 

The  computer  program  SHELL  was  used  exclusively  in  this 
research.  It  was  originally  developed  by  Dennis  (9)  for 
large  displacement/rotation  static  analysis  of  shell 
structures.  The  modified  Riks-Wempner  solution  algorithm 
was  added  to  the  program  by  Tsai  and  Palazotto  (26)  and  it 
was  extended  for  the  study  of  non-linear  vibrations  of 
cylindrical  shells  also  by  Tsai  and  Palazotto  (14).  The 
uniqueness  of  the  SHELL  program  is  that  parabolic  shear 
strains  are  assumed  to  vary  through  the  shell  thickness  and 
vanish  at  the  top  and  bottom  surface.  Although  classical 
plate  and  shell  theory  for  thin  isotropic  structures  ignore 
shear  stress  through  the  thickness,  it  is  not  appropriate  to 
do  so  with  composite  shells.  The  coupling  of  extensional, 
bending,  and  shear  strain  must  be  taken  into  consideration. 

It  also  should  be  noted  that  the  SHELL  program 
incorporates  material  linearity.  This  is  appropriate 
because  in  structural  applications,  materials  are  normally 
restricted  to  the  linear  elastic  region.  This  implies  that 
any  fiber  or  ply  breakage  during  large  rotations  and 
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displacements  is  ignored. 

SHELL  offers  two  strain  displacement  formulation 
options  for  a  cylindrical  shell  finite  element.  When  Dennis 
originally  wrote  the  code,  he  incorporated  the  Donnell  shell 
equations  along  with  the  non-linear  equations  he  developed 
in  order  for  comparisons  to  be  made.  Only  Dennis'  non¬ 
linear  option  is  used  in  this  research. 

The  analytical  approach  taken  in  this  work  is  to  use 
both  the  static  and  dynamic  analytical  capabilities  of  SHELL 
to  do  collapse  analysis  of  circular  shells  and  arches 
subjected  to  transverse  loading.  Comparisons  to  other  work 
is  done  when  applicable.  In  the  interest  of  semantics,  it 
should  be  noted  that  the  terms  buckling,  collapse,  and  snap- 
through  are  all  used  interchangeably. 

Previous  Work 

The  engineering  analysis  of  general  shells  goes  back  no 
further  than  the  start  of  the  twentieth  century.  Simmons 
(24)  notes  that  in  1920,  A.  E.  Love  published  a  set  of 
equations  for  the  midsurface  displacement  of  circular, 
cylindrical  thin  shells.  Perhaps  the  first  useful 
cylindrical  shell  equations  were  presented  by  Donnell  (10) 
in  1933.  He  determined  (both  analytically  and 
experimentally)  that  circumferentially  trigonometric 
deformations  with  small  wavelengths  allow  one  to  discard 
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several  terms  in  the  curvature  and  twist  equations.  The 
complexity  of  the  resulting  expressions  are  on  the  order  of 
the  Von  Karman  plate  equations.  The  primary  restriction  is 
that  the  ratio  of  shell  thickness  must  be  small  and  if 
rotations  of  the  midsurface  greater  than  fifteen  degrees 
occur  then  accuracy  suffers. 

Another  major  milestone  in  shell  theory  occurred  in 
1959  by  Sanders  (22)  .  He  assumed  that  transverse  shear 
strains  are  negligible  and  solved  for  the  corresponding 
transverse  shear  stress  resultants,  thus  partially 
incorporating  transverse  effects.  The  significant 
improvement  over  Love's  approximations  is  that  by 
considering  rotation  of  the  shell  normal,  but  neglecting 
rotations  about  this  normal,  Sanders'  relations  allow  for 
small  strain  free  rigid  body  motion. 

Since  finite  element  analysis  came  into  being  for  use 
in  structural  analysis,  many  researchers  have  developed 
finite  elements  for  cylindrical  shells.  Most  applicable  to 
the  present  research  is  the  work  done  by  Saber  and  Lock 
(21) ,  which  employed  a  solution  algorithm  capable  of  tracing 
the  post-collapse  behavior  by  incrementing  either  load  or 
displacement  to  avoid  numerical  singularity  at  critical 
points.  A  20  degree  of  freedom  curved  rectangular  element 
for  an  isotropic  material  was  used.  Most  theoretical  work 
done  up  until  this  point  in  time  involved  isotropic 
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materials.  With  the  advent  of  composites,  it  was  necessary 
to  extend  first-order  theory  to  laminated  materials.  Chung 
and  Widera  (4)  expanded  Donnell's  equations  for  use  with 
laminated  composites.  It  should  be  noted  that  the  limits  of 
the  applicability  of  the  theory  as  far  as  magnitude  of 
deformations  are  concerned  remained  the  same. 

All  work  reported  up  to  this  point  is  considered  first- 
order  theory.  The  development  of  a  higher  order  theory 
incorporating  transverse  effects  was  driven  by  the  need  to 
accurately  analyze  laminated  composite  shells  which 
typically  have  more  severe  transverse  stresses.  Dennis  (8) 
and  Reddy  and  Liu  (17)  have  both  presented  comparable  higher 
order  theories  incorporating  transverse  effects,  of  which 
Dennis'  is  the  basis  for  the  present  work.  This  theory 
allows  for  fully  non-linear  in-plane  strains,  but  only 
linear  transverso  shear  strain-displacement  relations.  The 
acceptability  of  a  linear  transverse  strain  field  in  a  non¬ 
linear  theory  is  justified  by  noting  that  transverse  effects 
are  small  compared  to  in-plane  effects,  thus  smaller  higher- 
order  terms  of  the  already  small  transverse  terms  are 
negligible  (9) . 

The  solution  method  for  geometrically  non-linear 
problems  by  the  finite  element  method  requires  recalculation 
of  the  structure's  stiffness  matrix  as  it  changes  during 
deformation.  The  Newton-Raphson  method  for  iteration  is 
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widely  used  in  such  computations.  However,  once  the 
structure  has  reached  its  peak  load,  it  has  essentially  zero 
stiffness  and  inverting  the  stiffness  matrix  as  required  for 
load-controlled  Newton-Raphson  iteration  to  converge  to  the 
equilibrium  solution  is  impossible.  This  problem  has  been 
circumvented  by  reformulating  the  equilibrium  equations  to 
step  by  displacement  rather  than  load  when  a  peak  load  was 
encountered,  and  using  the  usual  load  incrementing  method  to 
step  through  areas  of  zero  or  reversing  incremental 
displacement.  A  solution  algorithm  without  the  need  to  swap 
between  methods  was  developed  by  Riks  (18,19)  and  Wemper 
(27)  and  applied  to  structural  problems  by  Crisfield  (7) . 

In  it,  a  selected  arc  length  of  the  equilibrium  path  is 
incremented  rather  than  load  or  displacement.  Since  step 
size  does  not  decrease  to  zero  due  to  convergence  to  a  peak 
load  or  peak  displacement  point  as  the  former  methods  do, 
this  technique  steps  past  such  singular  areas  and  allows 
uninterrupted  tracing  of  the  equilibrium  path  (15).  Tsai 
and  Palazotto  (14)  have  applied  this  solution  algorithm  to 
the  SHELL  computer  code  used  in  this  research. 

Although  a  great  deal  of  research  has  been  done  into 
the  static  analysis  of  shell  structures,  the  development  of 
the  finite  element  method  for  nonlinear  shell  dynamic 
analysis  has  not  been  as  rapid.  Clough  and  Wilson  (5)  have 
done  work  into  nonlinear  vibrations  using  flat  plate 
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elements.  The  higher  order  transverse  shear  deformation 
theory  used  in  this  research  has  been  extended  for  the  study 
of  nonlinear  vibration  of  cylindrical  shells  by  Tsai  and 
Palazotto  (26) . 

Current  Work 

The  subject  of  this  research  is  the  dynamic  response  to 
a  transverse  point  load  of  a  laminated  cylindrical  shell  and 
arch.  The  static  and  dynamic  analysis  capabilities  of  the 
SHELL  computer  code  are  further  explored.  Both  dynamic  and 
static  analysis  of  the  post  buckling  response  of  cylindrical 
shells  are  compared.  No  experiment  or  comparable  numerical 
study  could  be  found  with  which  to  compare  results;  hence 
trends  and  general  conclusions  will  be  compared  with 
conclusions  in  some  of  the  studies  referenced  above.  The 
computers  used  were  an  Elxsi  6400  using  the  Unix  operating 
system  located  at  the  Air  Force  Institute  of  Technology  and 
a  Cray  Y-MP8/864  running  Unicos  located  at  the  Ohio 
Supercomputer  Center,  the  Ohio  State  University. 
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II.  THEORY 


In  order  for  the  reader  to  better  understand  the 
theoretical  basis  for  the  data  presented  later,  the  finite 
element  formulations  governing  the  SHELL  computer  code  is 
now  discussed.  The  geometry  and  assumptions  used  in  the 
model  will  be  presented  along  with  the  constitutive 
development,  the  strain-displacement  relations  and  the 
potential  energy  theory  used  to  develop  the  equations  of 
motion.  The  resulting  36  DOF  element  will  be  presented. 
Finally,  the  extension  of  the  code  to  dynamic  analysis  by 
Tsai  and  Palazotto  (14)  will  be  discussed. 

Bathe  and  Ho  (2)  stated  the  following  set  of  criteria 
for  a  desirable  shell  element: 

1. )  No  spurious  zero-energy  modes  should  exist,  so 

that  reliable  results  can  always  be  expected.  No 
numerical  fudge  factors  should  be  necessary, 
either. 

2. )  The  element  should  be  applicable  to  general  shell 

structures,  including  those  with  beam  stiffeners, 
cutouts,  intersections,  etc. 

3. )  The  element  should  be  cost-effective  for  linear  as 

well  as  nonlinear  static  and  dynamic  analysis. 

This  implies  that  the  degrees  of  freedom  is  held 
to  a  minimum.  It  should  allow  analysis  of  large 
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displacement  and  large  rotation  problems,  and 
materially  nonlinear  situations. 

Criterion  1  is  achieved  in  this  formulation.  The  code  has 
not  been  developed  beyond  the  research  stage,  so  application 
to  other  than  plate  and  cylindrical  shell  cases  modeled  with 
rectangular  elements  is  not  yet  possible.  However,  the 
underlying  theory  by  which  transverse  shear  is  incorporated 
is  not  restricted  geometrically,  so  criterion  2  is  somewhat 
met.  Material  linearity  is  assumed  in  the  formulation  of 
the  element,  so  criterion  3  is  only  partially  fulfilled. 

The  incorporation  of  through  the  thickness  shear  effects  in 
a  shell  structure,  while  maintaining  a  two  dimensional 
analysis,  satisfies  the  first  part  of  this  criterion. 

Geometry  and  Assumptions 

The  curvilinear  orthogonal  coordinate  system  and 
nomenclature  used  in  this  formulation  of  the  laminated 
cylindrical  shell  is  shown  in  Fig  (2-1) .  The  x-axis  lies 
along  the  straight  dimension  of  the  panel;  the  s-axis 
follows  the  circumference,  and  the  z-axis  is  everywhere 
normal  to  the  shell  middle  surface,  positive  toward  the 
center  of  curvature.  The  surface  formed  by  the  x  and  s  axes 
lies  in  the  center  of  the  thickness  of  the  panel,  so  the 
thickness  coordinate  is  negative  on  the  outer  surface  and 
positive  on  the  inner  surface.  Displacements  along  the  x,s 
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Figure  2-1.  Cylindrical  Shell  Geometry 


and  z  axes  are  u,  v  and  w  respectively.  In  the  early 
vectorial  development,  the  coordinate  0  is  used  instead  of  s 
for  the  purpose  of  generality  and  s  is  used  when 
specializing  to  the  cylindrical  geometry.  Since  the 
structure  analyzed  here  is  an  open  shell,  the  angle  0  is 
also  useful  for  describing  shallowness.  The  angle  <p 
specifies  the  orientation  angle  of  each  ply  in  the  laminate. 
Subscripts  denoting  stress  and  strain  orientation  are 
explained  in  Table  2-1  and  Fig  (2-2) . 
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Table  2-1.  SHELL  Contracted  Notation 


STRESS 


STRAIN 


EXPLICIT  CONTRACTED  EXPLICIT 


CONTRACTED 


CYLINDRICAL 


COORDINATES 


Figure  2-2.  Fiber  Reinforced  Lamina 
Definitions 
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Constitutive  Development 


In  this  section  the  stress-strain  relations  for  a 
laminate  of  arbitrarily  oriented  transversely  isotropic 
plies  will  be  developed.  The  main  difference  from  isotropic 
relations  is  the  summation  over  the  lamina  thickness  of  the 
directional  constitutive  equations  for  each  ply,  arriving  at 
the  total  laminate's  effective  stress-strain  relation.  This 
section  follows  the  references  by  Silva  (23)  and  Dennis  (8) . 

For  an  isotropic  material,  stress  a  and  strain  e  are 
related  as 

a  =  E  e  (2-1) 

where  E  is  Young's  modulus.  When  considering  the  general 
case  of  an  anisotropic  material,  Young's  modulus  can  differ 
with  different  load  orientations  so  this  equation  expands  to 

*-11  *-12  *-13  *-14  *-15  *-16 
*-21  *-*22  *-23  *-24  *-25  *-26 

*-31  *-32  *-33  *-34  *-35  *-36 

*-41  *-42  *-43  *-44  *-*45  *"46 

*"51  *-52  *-53  *-54  *-55  *-56 
*-61  *-62  *-63  *-64  *-6  5  *-66. 

where  C,..  is  the  stiffness  matrix  of  36  terms  defining  the 
stress-strain  relationship  for  loading  in  the  ith  direction. 
If  one  only  considers  the  energy  conserving  elastic  region, 
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the  matrix  C,- j  is  symmetric;  i.e.  C12  =  C21  and  so  forth, 
resulting  in  21  independent  terms.  In  the  case  of  a 
directional  fiber-reinforced  composite  material  (the 
material  in  this  research) ,  the  three  mutually  orthogonal 
planes  of  symmetry  decouple  shear  strains  from  normal 
stresses  and  vice  versa.  This  defines  an  orthotropic 
material,  which  has  only  9  independent  stiffnesses. 


°1 
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^13 
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0  ' 
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°2 

^12 

^22 

^"23 

0 

0 

0 

e2 

°3 

>  =  - 

^13 

^23 

*-33 
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0 

0 

►  i 

E3 

°4 

0 

0 

0 

*"44 

0 

0 

e. 

°5 

0 

0 

0 

0 

^55 

0 

e5 

°6 

0 

0 

0 

0 

0 

*-66. 

E6. 

(2-3) 


Also,  since  such  a  material  responds  equally  to  any 
direction  of  load  in  the  plane  perpendicular  to  the  fiber 
longitudinal  axis  (2-3  plane) ,  the  2  and  3  subscripts  are 
interchangeable.  This  behavior  is  called  transverse 
isotropy  and  further  reduces  the  number  of  independent 
stiffness  terms  to  seven.  Thus,  in  terms  of  engineering 
constants  E  (Young's  modulus)  and  v  (Poisson's  ratio)  the 
following  are  the  stiffness  terms; 
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i-vL 


cn~Ei  A 


E22  *-33  E2 


1  V12V21 


1  +V 


<'12  ^13  ^1V21  A 


23 


^2  3  E2 


V23+V12V21 


*'44=<-*23 
^S5=^31 
E66  =<^12 


(2-4) 


where  A  =  1— 2 v  12v2i~ v232- 2 v i2v2i v23 '  Ei  Youn<?'s  modulus  for 
loads  along  the  1  axis,  G12  is  the  shear  modulus  in  the  1-2 
plane,  and  v12  denotes  the  ratio  of  strains  e2/e 1  for  stress 
applied  in  the  1  direction. 

Due  to  the  fact  that  in  most  composite  materials  the 
plies  are  thin,  the  assumptions  of  plane  stress  (a3=a4=a5=0) 
is  usually  made  at  this  point.  However,  in  the  SHELL  code, 
non-zero  through  the  thickness  shear  stress  is  allowed,  thus 
a  "modified"  state  of  plane  stress  is  assumed  in  which  only 
<7j=0.  Solving  for  e3  in  Eq  (2-3)  after  applying  this 
assumption  yields 
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Applying  the  modified  plane  stress  assumption  to  Eq 
(2-3) ,  using  the  relation  of  Eq  (2-5) ,  and  rearranging  terms 
produces  the  lamina  constitutive  relation 


°1 

‘Oil 

012 

0 

0 

0 

o' 

El 

°2 

012 

022 

0 

0 

0 

0 

>  = 

0 

0 

066 

0 

0 

0 

1 

e6 

°4 

0 

0 

0 

044 

0 

0 

e4 

°5 

0 

0 

0 

0 

055 

0 

e5 

(2-6) 


where  the  Q^.  are  reduced  stiffness  coefficients  related  to 
the  C ' s  by 


Onsets 


G33 


(2-7) 


In  terms  of  engineering  coefficients, 


022  =  — 

CO 

0  ^21^ 

Q 

<0 

(2-8) 

012 

044=023 

055  =  G13 

where  o  =  1-  v12  v21. 

Finally,  in  order  to  analyze  a  stack  of  plies,  they 
must  all  be  referenced  to  a  global  axis  system  and  their 
effects  summed: 

(£>„],  inTfeA  <2~9> 
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where 


C2 

s 2 

-2  CS 

OlX 

012 

0 

[T]  = 

s2 

c2 

2  CS 

for 

0X2 

O22 

0 

cs 

-cs 

c2-s2 

0 

0 

£*66. 

[T] 


c  -s 
s  c 


for 


Qti 
0  fi 


0 

55 


and  c=cos($) ,  s=sin(*) .  With  this  transformation  the 
constitutive  relations  are 

(°A=  feA  <2-10> 

with  the  transformed  reduced  stiffnesses 

Qn  =  Q^cos4*  +  2  (Q12  +  20^)  sin2*cos2*  +  Q22sin4* 

Q12  =  (Qn  +  Q22  “  ^Q^)  sin2*cos2*  +  Q12(sin4*  +  cos4*) 

Q22  =  Qnsin4*  +  2  (Q12  +  2(2^)  sin2*cos2*  +  Q22cos4* 

Qi6  =  (Qn  ~  Q1Z  ~  2Q66)sin*cos3*  +  (Q12  -  Q22  +  2(2^)  sin3*cos* 

Q26  =  (Qn  "  Q\z  ~  20^)  sin3*cos#  +  (Q12  -  Q22  +  2Q,*)  sin*cos3* 

^66  =  (Q11  +  Q22  ~2Q12  “  2Q66)  sin2*cosz*  +  Q^sin4*  +  cos4*) 

Q44  =  Q44COS2*  +  Q55sin2* 

Q45  =  (Q44  “  Q55)  cos*sin* 

Q55  =  Q55COS2*  +  Q^sin2* 
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Strain  Displacement  Relations 


The  strain  versus  displacement  relations  in  the  SHELL 
code  take  into  account  the  geometric  nonlinearity  arising 
from  the  panel’s  curvature.  It  is  also  within  these 
relations  that  through  the  thickness  shear  effects  are 
incorporated  into  the  analysis.  This  development  follows 
Dennis  (8)  who  authored  the  original  code. 

The  incorporation  of  transverse  shear  effects  is 
accomplished  by  assuming  a  modified  state  of  plane  stress 
for  the  laminar  (previously  discussed)  in  which  a3=0  (and 
hence  e3=0)  but  u4  and  a5  are  allowed  to  be  small  non-zero 
values.  These  transverse  shear  stresses  are  assumed  to 
equal  zero  on  the  top  and  bottom  surfaces  of  the  shell,  and 
the  associated  strains  will  vary  parabolically  through  the 
shell  thickness. 

To  start,  one  looks  at  the  fully  non-linear  strain- 
displacement  relations  for  an  orthogonal  curvilinear 
coordinate  system  taken  from  Saada  (20)  and  presented  below. 


0UX 


.  hlU2  9hl 

h2  dy2 


+  hlU3  9hl 

h2  dy3 


+  1  (  dui  +  u2  dbi  +  u3  2 
2  dy1  h2  dy2  h3  dy3 


+  1  (  dU2  _  U1  dhl  )  2  t  1  /  dU3  _  U1  dhl  ,  2 

2  dy1  h2  dy2  2  dyt  h3  dy3 


(2-11) 
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v  =h  du2  h2u2  dh2  hziii  dhz 

22  2  dy2  *3  dy2  *x  dy1 

^  1  ,  dU2  A  U3  dh2  ^  U1  ^2  >  2 
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(2-13) 


(2-14) 
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(2-15) 


(2-16) 


The  u{ ’s  are  the  displacements  in  the  1,2, and  3  directions 
respectively  and  the  y{ 's  are  the  generalized  coordinates  in 
the  1,2,  and  3  directions.  The  expressions  for  e ..  are 
obtained  by  dividing  the  y^.  by  the  product  of  the 
coordinate  system  scale  factors  h;hj.  It  is  important  to 
keep  in  mind  the  contracted  notation  discussed  in  Table  2-1. 
As  stated  earlier  e3  is  assumed  to  be  zero  and  only  the 
linear  terms  are  kept  for  the  transverse  shear  strains  e4 
and  e5.  Thus  we  see  that  (8) : 
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(2-17) 


1 

e4  =  —  [u3 , 2  +  h2u3,3  -  u2h2,3l 
es  =  -j-  [u3/1  +  h^,  3  -  uxi3lf3] 

■"i 

where  the  coordinate  system  scale  factors  for  the 
cylindrical  geometry  used  in  this  work  are  h^l,  h2=l-  z/R, 
and  h3=l. 

The  displacement  equations  in  the  thickness  variable  z 
which  permit  the  incorporation  of  the  desired  through  the 
thickness  features  are  (8) : 


u(x, 0, z ) =  u° 
v(x,Q,  z)  =  v° 
w(x, 0) =  w 


+  zty1  +  z24>x  +  z3yx  +  z48x 


1 


i? 


+z\lf2+z2<t>2+z3Y2+z492 


(2-18) 


where  u°,  v°,  w,  i|r. ,  <p{,  y,»  and  0,-  are  functions  of  the 
coordinates  x  and  0.  The  displacements  u°  and  v°  are  of  the 
shell  middle  surface.  The  transverse  displacement  w  is  the 
same  throughout  the  thickness  since  transverse  normal  strain 
is  assumed  negligible.  The  terms  are  rotations  of  the 
surface  normals  in  the  x  and  s  planes,  and  <p- ,  and  ©j 
are  to  be  found  by  applying  the  assumptions  that  transverse 
shear  stresses  a4  and  a5  are  zero  in  the  shell  surfaces. 
Substituting  the  equations  for  v  and  w  into  Eq  (2-17)  for  e4 
yields  (8) : 
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~~R  +  * 2  +  2z*a  +  3z2y2+  4z3e2j 

(2-19) 

+  l(v°(l-|)  +  z*|f2  +  z24>2  +  z3y2  +  z492 

For  zero  transverse  shear  stress  at  the  surfaces,  the 
associated  strain  will  also  be  zero.  Enforcing  this 
condition  by  substituting  ±h/2  for  z  in  Eq  (2-19)  and 
setting  both  resulting  expressions  to  zero  and  hence  equal 
to  each  other,  one  obtains  by  then  solving  for  the  unknown 
variables  the  following  (8) : 


4>2  =  0 


«2 


Y2 

2  R 


(2-20) 


Note  that  an  h/R  value  of  1/5  (quite  high  for  practical 
aerospace  shell  geometries)  allows  for  the  neglection  of  the 
h/R  term  in  the  left  side  of  Eq  (2-20) .  Replacing  02,  02 
and  y2  in  Eq  (2-19)  with  the  results  in  Eq  (2-20)  one  can 
find  the  transverse  shear  strain  in  terms  of  transverse 
displacement  w  and  rotation  y2  (8)  : 


64  1-z/R 


( w 


*4-  4P 


+  8 


3  h2R 


(2-21) 
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It  should  be  noted  that  there  are  1/R  terms  to  be  neglected. 
Therefore,  the  transverse  shear  strain-displacement 
relations,  c4,  along  with  a  similar  analysis  with  e5 
(simpler  since  h,  =  1)  ,  become  (8)  : 


«4  = 


1 

1-z/R 


(W'2 


+ 


(2-22a) 


e5= 


(w, 


+ 


(2-22b) 


By  replacing  <p2,  02  and  y2  in  Eg  (2-18)  it  is  found  that  the 
displacement  equations  are  (8) : 


u(x,e,z)  =  u°  +  z^1  -  +  w>i) 


v{x,  0,  z)  =  v° 
w(x,  0)  =  w 


1_l]+  zfz"  ^z3(fz+  w’2) 


(2-23) 


At  this  point,  one  can  note  that  this  formulation  provides 
seven  degrees  of  freedom:  u;  v;  w;  wM;  w,2;  \(»1  and  if2. 

Now  that  the  displacement  equations  which  incorporate  a 
parabolic  through  the  thickness  shear  stress  distribution 
have  been  developed,  the  in-plane  kinematic  equations  are 
derived  for  the  shell  middle  surface.  The  fully  general 
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strain  displacement  relations  are  quite  extensive;  the 
theoretical  development  incorporating  them  can  be  found  in 
Dennis'  large  displacement  -  moderate  rotation  general  shell 
development  (8) .  The  general  strain  displacement  relations 
of  Eq  (2-11)  and  (2-16)  with  the  kinematics  of  Eq  (2-23) 
will  give  the  inplane  shell  strain  displacement  relations. 
These  expressions  can  be  specialized  for  a  shell  geometry  of 
interest  by  defining  the  scale  factors,  h^  using  Eq  (2-24) 
(8)  . 

hj  =  aj  (1  -  z/R2) 

h2  =  CC2(1  -  z/R2)  (2-24) 

h3  =  1 

where  ay2  =  gYY  (no  sum)  ,  which  are  called  the  metric 
coefficients  for  the  orthogonal  curvilinear  coordinate 
system. 

By  substituting  Eq  (2-24)  into  Eq  (2-11)  -  (2-16)  and 
carrying  out  the  indicated  differentiation,  one  obtains  for 
the  inplane  strain-displacement  equations  (8) : 
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(2-25a) 


ei  =  =  e°  +  ziX  +  *2Xi 

hi 

+  z3Xi  +  z4xJ  +  z6xl 


=  111  =  e°  + 


*X2  +  ^2X2 


+  z3xl+  z*xi+  z6xl 


(2-25b) 


e£  =  4^  =  fig  +  zxe  +  z2xl 


h1h2 


+  *3xl  +  24X6  +  z6Xs 


(2-25C) 


where  the  ej°  and  the  Xj1  terms  (j=l,2,6;  i-1,2,3,4,6)  are 
functions  of  the  displacements  and  the  scale  factors,  and 
can  be  found  in  Appendix  A  of  (8) .  It  is  noted  that  the 
superscripts  on  the  x j 1  are  not  exponents.  They  are 
individual  strain  components  that  correspond  to  the  power  of 
z  that  multiplies  it,  and  the  subscripts  on  Xj'  indicate  the 
strain,  e,,  e2,  or  e6  that  these  components  correspond  with. 

The  following  equivalent  representation  of  the  strains 
is  conducive  to  the  matrix  operations  of  the  potential 
energy  formulation  in  the  next  section.  The  inplane  strain 
displacement  relations  are  represented  by  (8) : 
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(2-26a) 


=  e°  +  zpXip  ;  (P= 1/  •  •  #7) 

where 

€?  =  U, 2  +  -|  (u,i  +  V,\  +  W,2) 


e2  =  €2  +  zpX2p  ;  (P=l/  •  •  /  7) 

where 


e°=v/ 


2-H(u-+  v-+ 


/  2 


*i) 

i?2  i?2 


+  VW'  2  ^  V'2W 
R  "  R 


(2-26b) 


e6  =  e6  +  ^px6P  ;  (p=i/  •  •  /  7) 

where 


„o  _ 


e6  =  u,2  +  v,  x+  u#1u,2+  vflv,2+  wflw,2 
+  4  (VW,  2  -  v,  xw) 


(2-26c) 


Note:  The  higher  order  bending  strain-displacement 

relations,  %\p  (i=l,2,6;  p=l,..,7),  are  completely  shown  in 
(8)  . 

The  transverse  shear  strains  (in  similar  form)  are  (8) 

^4  =  w,  2  +  i|f2 

e4=  zpxip  ;  X42  *  3 k(w,2  +  i|f2)  (2-27a) 

X4p(p=l,3,4,5,6,7)  =0 
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e5  =  v/:l  +  i|r1 

e5=  e5+  zpx5p  ;  X52  =  ^k(wtl  +  i|r1)  (2-27b) 

X5p  (P=l ,  3 , 4 , 5 , 6 , 7 )  =0 


This  allows  assembly  of  the  strain  displacement  equations 
into  a  matrix  format: 


*1 

e3 


ro 

r1 

Xu 

X12 

Xl3 

Xl4 

X15 

Xl6 

Xl7 

r° 

■  + 

X21 

X22 

X23 

X24 

X25 

X26 

X27 

Lo 

r3. 

X61 

X62 

X63 

Xe4 

X65 

X66 

X67 

+ 


X42 

X52 


(2-28) 


In  a  general  expression, 

{%}  =  {^°}  +  [X]  (Z) 


(2-29) 


It  is  worth  noting  that  these  kinematics  avoid  the 
common  pitfall  of  shear  locking,  wherein  the  model  becomes 
artificially  stiff  as  the  shell  thickness  is  decreased. 
This  is  a  problem  with  finite  element  formulations  which 
incorporate  constant  or  linearly  distributed  shear  strain 
through  the  thickness,  and  it  necessitates  the  use  of  a 
correction  factor.  However,  examination  of  the 
compatibility  relations  associated  with  the  strain 
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displacement  equations  developed  in  this  section  shows  that 
the  terms  associated  with  transverse  shear  drop  out  as 
thickness  is  reduced  to  zero  (8) . 

Note  the  kinematics  developed  in  this  section  are 
specifically  for  the  cylindrical  geometry.  The  analysis 
accomplished  for  this  work  uses  these  large  displacement  - 
moderate  rotation  formulations. 

Potential  Energy 

The  shell  potential  energy  is  the  sum  of  the  internal 
strain  energy  and  the  work  done  by  external  forces: 

J|p  =  E  +W  (2-30) 


where  the  internal  strain  energy  is  given  by 

E  =  [  [  dz  dQ  (2-31) 

JQJ  h  2 


where  n  represents  the  shell  middle  surface.  The  internal 
strain  energy  E  is  composed  of  in-plane  and  normal  terms 
(set  E1 )  and  transverse  shear  terms  (set  E2)  .  Expanding  the 
expression  for  E  by  inserting  Eqs  (2-27)  through  (2-29)  into 
Eq  (2-31)  gives 
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El  =  iUjQ“(e0l  +  ZPXl*>)2  +  Q22^2  +  ZPX2p)2 

+  2Q12(el+z*x iP)  (e2°+z%r) 

+  066  <e°  +  zPX6p)2  (2-32) 

+  2016(e?+zpXiP)  (eg  +  zrx6r) 

+  2p26(e°+zpX2P)  (Gs  +  zrx6r)]  dz  dQ 

Bi  =  if0fJ°“(e‘*z2x->>* 

+  p55(e^zpx52)  2  (2-33) 

+  2045(C4+z2x42)  (e°+z2x52)]  dz  dQ 

where  p, r=l , 2 , . . . , 7 .  Integrating  the  z  over  ±  h/2  yields 
the  equation  for  strain  energy  as  a  function  of  the  middle 
surface  only,  which  is  the  desired  formulation  for  this 
problem.  A  further  simplification  performed  in  the 
development  of  the  SHELL  program  is  one  of  symmetry  in  ply 
layup.  This  results  in  the  cancellation  of  elasticity 
arrays  which  are  multiplied  by  odd  powers  of  the  transverse 
coordinate  z.  Rearranging  yields  the  final  forms  of 

E1  =  -  f  {g°)T[A]{g°}  dQ  (2-34) 

2  J  Q 
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E2  =  If  (  {g0)T[A]{&0)  +  2  {g’°}r[D]  [X] 
2  2Jq 


(2-35) 


+  [X]T[F]  [X])  dd 

where  stiffness  matrices  { [A,  D, F] )  =  /h  [Q]  {1,  z2,  z4}  dz. 
The  36  DOF  Element 

The  SHELL  computer  code  uses  a  36  DOF  element 
developed  by  Dennis  (8)  and  \  tured  in  Fig  (2-3) .  The 
seven  degrees  of  freedom  at  each  corner  node  were  found  in 
the  previous  section  to  be  u,  v,  w,  wM,  w,2,  y1  and  y2.  The 
mid-side  nodes  have  only  the  two  inplane  degrees  of  freedom 
u  and  v.  C°  continuity  is  required  of  all  but  w  and  its 
derivatives  which  require  C1  continuity.  Lagrangian 
bilinear  interpolation  is  used  for  u,  v,  y1  and  y2  and  non 
conforming  Hermitian  interpolation  is  used  for  w,  wM  and 
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Figure  2-3.  SHELL  36  Degree  of  Freedom  Element 


The  assumed  transverse  displacement  equation  for  the 
cylindrical  shell  element  is 

w(x,  s)  =a1  +  a2x  +  a3s  +  aAx2  +a5xs  +  a6s2  +  a7x3  +  a8x2s 

+  a9xs2  +  a1Qs3  +  anx3s  +  a12xs3  (2-36) 

This  rectangular  element  may  be  isoparametrically  scaled  and 
oriented  such  that  the  longitudinal  coordinate  x  -*  i ,  the 
circumferential  coordinate  s  -*  r) ,  and  the  side  lengths  are 
scaled  by 

Z  =  -  n  =  |  (2-37) 

a  o 
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where  a  and  b  are  the  scaled  element  half  dimensions  in  the 
5  and  rj  directions  respectively.  Now  the  element  transverse 
displacement  equation  can  be  written  as: 


w{x,  s) 


lXi  X2  X3  Xj 


(2-38) 


where  (q}kT=  {w  wM  w,2)  for  the  kth  corner  node,  and  the  x,'s 
are  Hermitian  shape  functions  (6): 


xh 


Xki 

xk2 

xk3 


(2+5*5 +ti*Ti-s2-n2) 

-§«,(i*{,5>2«,5-i>  (i*<i*<i> 

-|i|,<i<-{,?>  (<1,11-1)  (i*n,<i)2 


(2-39) 


The  formulations  for  DOF's  u,  v  and  i|rf  is  a  simple 
Lagrangian  form,  and  the  midside  nodes  5-8  are  also 
included: 


0  0  0 
Q\  0  0 

0  Nx  0 
0  0  Nx 


qa  o  o  o  os 

0  £>4  0  0  0 

0  0  Nt  0  0 

0  0  0  Af4  0 


0 

0 

0 


QTi 

^2 

^3 

Qs’ 

<?6 

<3S 
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where  (q)kT=(u  v  ^  ^r2 >  for  nodes  1  through  4  and  {u  v}  for 
nodes  5  through  8,  and  N  and  Q  are  linear  and  quadratic 
(respectively)  Lagrangian  shape  functions  (6) : 

Nk  =  -Jd+W)  a+Ti*n) 

£>*=-| (1+?*?)  (l+r\kt\)  (ZkZ+T)kT)-l)  , 

k=l ,2,3,4  (2-40) 

Qk=^  (1~V)  (l+n*!))  *  k=6,8 
Qk=\  (1-n2)  (1+^0  ,  k=5 , 1 


The  overall  equation  defining  the  isoparametric 
discretization  of  continuum  displacements  into  nodal 
displacements  is  obtained  by  merging  the  above 
relationships : 


|u|  =  Civ]  {g) 

(7X1)  {1x36)  (36 Xl) 


(2-41) 


Static  Finite  Element  Solution 

The  solution  to  the  static  finite  element  problem 
involves  finding  the  equilibrium  state  between  applied  load 
and  structural  response.  This  state  can  be  determined  by 
finding  where  the  variation  of  the  system  potential  energy 
is  zero.  Recall  Eq  (2-30) : 
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Ilp  ~E*U 


(2-30) 


where  now  internal  strain  energy  can  be  represented  by 

E  =  ±qt  q  =  —  qT  K  q  (2-42) 

3  [  2  3  2 

where  q  is  a  column  array  of  nodal  displacements,  K 
includes  the  constant  stiffness  terms,  N1  includes  the 
stiffness  terms  linear  in  displacement  and  N2  includes  the 
terms  quadratic  in  displacement.  The  external  work  can  be 
represented  as 

w  =  -qT  X  P  (2-43) 

where  P  is  a  column  vector  of  applied  nodal  loads  and  X  is  a 
multiplier  which  will  be  discussed  in  the  next  section  on 
the  static  solution  algorithm.  Equilibrium  is  defined  as 
the  state  at  which  internal  strain  energy  and  external  work 
balance,  thus  potential  energy  will  be  at  a  relative 
minimum.  This  point  can  be  found  by  substituting  Eqs  (2-42) 
and  (2-43)  into  Eq  (2-30)  and  taking  the  first  variation: 

8IIp  =  bqT  [K  q  -  Xp]  =  0  (2-44) 

Since  displacements  <5q  are  nonzero  for  all  but  the  trivial 
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solution,  the  bracketed  expression  must  be  zero  for 
equilibrium.  It  is  a  function  (called  f)  of  q  and  k: 

K  q  -  kp  =  0  =f(q,k)  (2-45) 

Since  K  varies  with  load  and  displacement,  a  numerical 
iterative  solution  is  used  to  solve  Eq  (2-45)  incorporating 
the  Newton-Raphson  method. 

Static  Solution  Algorithm 

For  static  analysis,  the  SHELL  code  uses  the  solution 
technique  advanced  by  Riks  (18,19)  and  Wempner  (27)  and 
demonstrated  by  Crisfield  (7)  of  incrementing  a  desired  arc 
length  along  the  load-displacement  curve  while  solving 
iteratively  via  the  Newton-Raphson  method.  Called  the 
modified  Riks -Wempner  method,  and  added  to  the  SHELL  code  by 
Tsai  and  Palazotto  (26) ,  this  allows  for  tracing  of  the  load 
displacement  response  through  both  load  reversing  (snap 
through)  and  displacement  reversing  (snap  back)  critical 
points,  so  the  most  complex  behavior  can  be  continuously 
followed. 

The  essence  of  the  Riks-Wempner  method  is  that 
neither  load  P  nor  displacement  q  is  independently 
controlled;  rather  a  selected  "arc"  length  As  (actually  the 
chord)  of  the  load-displacement  curve  is  incremented.  The 
equilibrium  condition  is  found  which  satisfies  the  relation: 
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Agi+1  •  Agi+1  +  AX*+1  P-P  =  As2 


(2-46) 


where  A|>1  is  the  incremental  displacement  for  step  i+1,  and 
Ali+1  is  the  fraction  of  load  P  applied  at  step  i+1.  The 
effect  of  the  constraint  equation  (2-46)  is  that  each 
subsequent  step  solution  is  searched  for  on  an  arc  of  radius 
As  from  the  current  solution.  The  initial  value  for  the 
quantity  AX  is  specified  with  the  problem  input  data. 

To  apply  the  Newton-Raphson  method,  the  first 
variation  of  Eq  (2-45)  is  taken  and  applied  at  step  i: 

Kt  Ag,  =  6XSP  -  Jf(gi,Xi)  (2-47) 

where 

=  5<7ii  +  (2-47a) 

the  out-of-balance  term  is 

bqn  =  -Ktx  (2-47b) 

the  linear  term  is 

bqi2  =  Kt  P  (2-47c) 

and 
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Kt  =  [K  +  Nx  *  N2] 


(2-47d) 


The  following  briefly  describes  the  operation  of  the 
algorithm  on  load  step  n,  i  denotes  an  iteration  at  step  n 
along  the  solution  path,  and  Fig  2-4  provides  a  simplified 
illustration  of  the  algorithm. 

a. )  The  tangent  stiffness  matrix  KT  at  the  current 

deformed  geometry  is  determined. 

b. )  The  linear  incremental  displacement  S-2  is 

computed  (Eq  2 -4 7c) . 

c. )  The  first  iteration  computes  Aqt  =  AA.,6^,  with 

AA,1=AA.n.1,  or  a  user  defined  value  if  n=l  (first 
increment) .  The  parameter  AX  indicates  the 
fraction  of  total  load  to  be  applied  at  the  first 
increment . 

d. )  The  constraint  equation  (2-46)  is  solved  for  As. 

The  load  term  is  often  ignored  in  this  relation, 
since  the  load  and  displacement  values  typically 
differ  by  many  orders  of  magnitude,  which  can 
cause  numerical  difficulty.  Convergence  to  a 
solution  is  not  hindered  by  ignoring  load  at 
this  step. 

e. )  KT  is  updated  at  q  =  q^  +  Aqj . 

f. )  The  out-of-balance  displacement  5n  is  computed 

(Eq  2-47b) . 
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g. )  Eqs  (2-47a)  in  (2-46)  are  solved  for  SXjr  which 

will  have  two  roots  51^  and  5Ai2  due  to  the 
quadratic  feature. 

h. )  is  selected  from  6An  and  <5Aj2,  by  criteria 

detailed  in  Crisfield  (7),  to  ensure  that  the 
solution  path  does  not  return  to  equilibrium 
points  previously  obtained. 

i. )  Displacement  increment  Aqj+1  =  Aqj+5q).  and  load 

factor  increment  AA1+1  =  AAj+SAj  are  updated. 
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j.)  When  Aqs  and  AXj  computed  on  successive  iteration 
(steps  e-i)  differ  by  less  than  a  selected 
convergence  tolerance,  the  step  n  solution  has 
been  found: 


Agi  -  A qn 


(2-48) 


At  the  completion  of  each  step,  cumulative  displacement 
and  load  factors  are  computed: 


Qn  =  Qn-X  + 
"  ^n-1  + 


(2-49) 


This  technique  can  follow  an  equilibrium  path  which 
progresses  in  any  direction,  by  solving  for  negative 
incremental  displacement  or  load,  or  both.  This  allows  one 
to  automatically  solve  for  even  the  most  convoluted 
nonlinear  equilibrium  curve,  without  concern  for  the 
singularities  at  critical  load  or  displacement  points. 

The  efficiency  of  the  algorithm  is  improved  by  scaling 
the  step  (n+1) *  s  target  As  length  by  the  ratio  of  a  user 
selected  desired  number  of  iterations  to  the  number  of 
iterations  to  converge  to  step  n.  Thus,  in  near-linear 
parts  of  the  load  displacement  curve,  wide  spacing  of 
solution  points  is  allowed;  when  the  curve  rounds  corners 
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(around  limit  points)  the  method  has  to  iterate  more  to 
converge  and  hence  As  is  reduced  until  the  curve  straightens 
out  again. 

The  inability  to  solve  for  equilibrium  at  a  limit  load 
point  (as  occurs  in  other  techniques)  is  circumvented  due  to 
the  nature  of  the  stepping  method.  The  technique  shoots  a 
tangent  from  the  current  equilibrium  point,  then  searches  an 
arc  about  the  tip  for  the  next  solution  point,  so  the  exact 
limit  load  point  is  almost  always  skipped  over.  Auxiliary 
equations  can  be  programmed  to  enable  a  more  exact 
determination  of  the  critical  point  (18). 

Dynamic  Equations  of  Motion 

Tsai  and  Palazotto  (14)  have  extended  the  SHELL  code 
for  the  study  of  nonlinear  vibration  in  cylindrical  shells. 
The  equations  of  motion  for  the  cylindrical  shell  are 
derived  via  Hamilton's  principal  where  the  variation  in  the 
time  integral  of  total  energy  is  set  to  zero: 

6  [*'  (E-T-W.)  =  0  (2-50) 

where  E  is  the  potential  strain  energy,  T  is  the  kinetic 
energy  and  We  is  the  work  of  the  external  forces. 

The  variational  components  6E,  6T,  and  <SWe  for  a 
laminated  plate  or  shell  are  given  by: 
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6E=f  T.  fCk 

ioAlJ<~  (2-51) 

+  Pjk)  6  Ujk)}  cfCdQ 

and 

L 

6 T=[  £  [k  {p'^vj^Svj*'}  dCdQ  (2-52) 

•’  °  Jtk-i  1 


and 


bWe  =  [  Fib  Uj  dQ 

e  JQ  J  3 


(2-53) 


where  i,j=l,2,3,  Ck.}  and  Ck  are  positions  at  the  bottom  and 
top  surface  of  the  k-th  layer,  n  is  the  domain  of  the 
neutral  surface,  a?j<k),  6ejj.<k),  Pj<k>,  6Uj(k>,  Vj(k),  p(k),  c(k),  are 
stress  tensor,  virtual  strain  tensor,  body  force  vector, 
virtual  displacement  vector,  velocity  vector,  mass  density, 
damping  coefficient  for  the  k-th  layer,  and  Fj  is  the 
external  force  vector  respectively,  L  is  the  total  number  of 
layers  in  a  laminated  shell. 

From  Eq  (2-50) ,  the  finite  element  formulation  is 
derived: 


MUi2)  +  CUa)  +  K(U{0))Ui0)  =  Fit) 


(2-54) 
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where  K  is  the  stiffness  matrix  for  large 
displacement/rotation  of  cylindrical  shells  which  was 
derived  earlier  and  presented  again  in  slightly  different 
form: 


K(U(0)) 


KAU(0)) 
Kn  +  — - — r - 


K2(U(0) ) 
+  — - - 

3 


(2-55) 


where  U<0),  U(1),  U<2)  are  the  displacement,  velocity,  and 
acceleration  vectors  at  each  nodal  point,  Kg  is  a  constant 
stiffness  matrix,  K1  is  a  stiffness  matrix  related  to  linear 
displacement,  K2  is  a  stiffness  matrix  related  to  quadratic 
displacement,  and  P(t)  is  the  external  load  applied  at  each 
node.  The  consistent  mass  matrix  is  obtained  as  described 
in  reference  (14)  by  substituting  Eq  (2-41)  into  Eq  (2-52) 


M 


L 

=  f  T  [  p<*>  [N]t[R]t[R]  [N]  dCdQ  (2-56) 

Jq  k=l  l 


where  [N]  is  the  overall  matrix  of  shape  functions  described 
in  Eqs  (2-39)  and  (2-40)  and  [R]  is  a  matrix  containing 
coordinate  system  scale  factors  as  shown  below: 


[K] 


10  0  Jc(3  0  ieC3+C  0 

0  1-A  0  0  Jci3  0  K3  +  C 

R 

1  0  1  0  0  0  0 
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where  k  =  -4/3h2,  h  is  the  shell  thickness,  R  is  the  radius 
of  the  shell  and  C  is  the  position  along  the  z  axis  or  the 
surface  normal.  The  damping  matrix  is  obtained  in  a  similar 
manner. 

L 

C  =  f  T  [k  c(k)  [MT[R]T[RUrf)  (2-57) 

jQ  k^l  Jik-' 


Solution  Technique  for  Dynamic  Problems 

To  solve  the  nonlinear  dynamic  problem  SHELL  uses  the 
beta-m  method  which  is  a  generalization  of  Newmark's  time 
marching  integration  scheme.  It  provides  a  general  single 
step  algorithm  applicable  to  initial  value  problems  and  is 
specialized  by  specifying  the  method  order  m  along  with  m 
integration  parameters,  f3Q,P,  ,/32, . .  .  ^or  a  particular 

choice  of  m,  the  integration  parameters  provide  a  subfamily 
of  methods  which  control  accuracy  and  stability.  The  finite 
difference  scheme  is  not  required  in  this  method.  This 
provides  a  better  way  for  computer  programming  than  the 
regular  methods  such  as  Newmark's  method.  The  beta-m  method 
is  defined  by 

u£\  =  gk  +  (2-58) 
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where 


Qk 


m 

-h 


U{nj)hj-k 

(j-k)  ! 


(2-59) 


and 


*>k  - 


(m-k)  ! 


(2-60) 


and  k=0, 1, . . . ,m,  is  defined  to  be  equal  one,  and  h  is  the 
time  increment  for  each  time  step.  The  method  order  m 
implies  that  Un(m)  is  the  highest  derivative  to  be  retained 
(m=2  is  Newmatk) .  For  this  research  m=2  for  all  analysis. 
Also,  throughout  this  work,  0O-J.5,  ^=0.5  and  (3y=1.0  as 
suggested  by  Fntona  and  Zienkiewicz  (12)  for  an 
unconditionally  stable  analysis.  Substituting  Eq  (2-58) 
into  Eq  (2-54)  at  time  tn+1 ,  results  in 


[. b2M+  i»xC+  b0K(  qQ  +  A  U{m) )  ]  A  U(m)  = 

Pn> i  "  {MQ2+  Cq1+  K(q0+  b0MJin>) )  q0) 


(2-61) 


where  is  the  applied  load  at  time  trHl ,  b0,b.,,b2,  are 
scalers  dependent  on  the  integration  parameters  as  shown  in 
Eq  (2-60),  and  q0,q1,q2  are  history  vectors  known  at  time  tn 
as  shown  in  Eq  (2-59) .  Eq  (2-61)  results  in  a  set  of 
nonlinear  algebraic  equations.  The  Newton-Raphson  iterative 
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method  is  adopted  here  to  solve  Eq  (2-61) .  At  each  time 
step  the  following  assumption  is  made 


AujTl  =  A  ujm)  +  6uJm) 


(2-62) 


where  i  is  the  iteration  number.  Applying  the  Newton- 
Raphson  method  to  Eq  (2-61)  and  using  Eq  (2-62)  yields 

[b2M+  bxC+  b0KT(g0+  b0At/iJB)]  bui  = 

Pn+ 1-  M{q2+  b2Aulm)}~  C  {gx+  Jb1At/jJB)}  (2-63) 

-  k(q0+  bQAu\m) ){<70+ 


where  KT  is  the  updated  tangential  stiffness  matrix. 

Eq  (2-63)  is  then  solved  by  the  following  algorithm: 

1.)  Given  U(0),  U' c1) ,  .  .  . .  at  time  t  ,  the  solution 

at  time  trr+1  is  desired. 

2*)  q0, q1 ,...., q^,  are  calculated  from  Eq  (2-59). 

3. )  Given  AUj1"0  and  Un+1<m)  from  the  i-th  iteration,  the 

right-hand  side  of  Eq  (2-63)  is  obtained  along 
with  the  updated  tangential  stiffness  matrix. 

4. )  Eq  (2-63)  is  used  to  solve  for  the  primary 

unknown, 

5. )  The  update  of  the  solution  vector  AUi+1(m)  is 

calculated  from  Eq  (2-62) . 

6. )  Um1<0),  Un+1(1>, ....  ,UrH.1(m)  is  updated  for  the  i+l-th 

iteration  from  Eq  (2-58) . 
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7.  )  If 


i+i 


N 


s  e 


(2-64) 


then  the  algorithm  goes  back  to  step  (1)  for  the 
next  time  step,  otherwise,  it  goes  to  step  (3)  for 
the  next  iteration.  1  is  the  degree  of  freedom 
number  and  L  is  the  total  number  of  dof's  in  the 
finite  element  model,  e  is  the  given  convergence 
criterion. 

The  above  procedures  are  terminated  when  tn  is  equal  to  the 
given  time  duration. 
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III.  RESULTS /DISCUSS IQ 


This  section  presents  results  of  analysis  done  on 
various  problems  to  further  evaluate  the  nonlinear  dynamic 
response  of  a  composite  cylindrical  panel.  A  static  snap- 
through  analysis  was  accomplished  on  an  isotropic  arch  using 
the  Riks  technique  and  compared  to  an  analysis  done  by 
Dennis  (8)  using  a  displacement  control  method.  Also,  a 
static  analysis  of  an  arch  is  compared  to  work  by  Belytschko 
and  Glaum  (3) .  Then  the  dynamic  response  of  an  isotropic 
panel  is  compared  to  some  present  results  by  Clough  and 
Wilson  (5) .  Finally,  the  majority  of  work  presented  shows 
both  a  laminated  composite  arch  and  shell  subjected  to  a 
pre-buckling  and  a  post-buckling  load  and  the  time  dependent 
behavior  compared  to  a  static  analysis. 

It  should  be  noted  that  all  finite  element  analysis 
done  in  this  research  took  advantage  of  using  1/4  panel 
symmetry  in  order  to  reduce  the  number  of  elements  and 
refine  the  mesh  size.  The  acceptability  of  using  symmetry 
with  the  SHELL  code  is  discussed  in  reference  (23)  by  Silva. 

Simple  Supported  Isotropic  Arch 

A  static  snap-through  analysis  of  a  simple  supported 
circular  arch  was  accomplished  where  the  geometrical  and 
material  quantities  are  shown  below. 
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E  =  10.0E6  psi 


v  =  0 

R  =  100  in.  (radius) 
h  =  1.0  in.  (thickness) 
w  =  1.0  in.  (width) 

0  =  53.13  deg 

The  mesh  incorporated  46  identical  elements  modeling  1/4  of 
the  arch  with  an  aspect  ratio  of  approximately  four  to  one 
as  shown  in  Figs  (3-1)  and  (3-2). 
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Figure  3-2.  Mesh  for  Isotropic  Arch 


The  loading  is  applied  radially  inward  exactly  at  the 
center  of  the  arch.  This  results  in  a  symmetric  response 
where  the  center  of  the  arch  displaces  only  radially.  The 
solution  by  Dennis  (8)  was  obtained  by  incrementing 
components  of  displacement  and  Fig  3-3  shows  a  comparison  of 
the  center  load  displacement  of  the  arch. 
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Figure  3-3.  Force  vs  Displacement  of  Center  Node,  Simply 
Supported  Isotropic  Arch 


As  can  be  seen,  both  analyses  predict  the  same  pre¬ 
buckling  load-displacement  curve  but  the  Riks  algorithm 
allows  for  a  tracing  of  the  equilibrium  path  to  a  new  post- 
buckling  solution.  The  Riks  method  steps  past  singular 
areas  and  allows  uninterrupted  tracing  of  the  equilibrium 
path. 
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Clamped  Isotropic  Arch 


A  static  collapse  analysis  of  a  clamped  arch  was 
accomplished  and  compared  to  work  done  by  Belytschko  and 
Glaum  (3) .  Geometric  and  material  quantities  are: 

E  =  10.0E6  psi  h  =  0.1875  in 

v  =  0  w  =  1.0  in 

R  =  113.114  in.  0  =  17.22  deg 

The  solution  uses  the  fully  nonlinear  strain-displacement 
relations  of  the  SHELL  code.  Mesh  sizes  of  10  and  20 
elements  were  used  with  almost  identical  results  thus 
convergence  is  assumed.  Belytschko  and  Glaum  reach 
convergence  for  their  fully  modeled  arch  when  the  mesh  size 
is  decreased  from  16  to  32  elements. 


Figure  3-4.  Clamped  Isotropic  Arch 
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Figure  3-5.  Mesh  for  Clamped  Isotropic  Arch 

The  loading  is  again  applied  radially  inward  at  the  center 
of  the  arch  which  results  in  the  center  of  the  arch 
displacing  only  radially.  Fig  (3-6)  shows  a  comparison  of 
the  two  load-displacement  curves. 
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Figure  3-6.  Load  vs  Displacement  of  Center  Node, 
Clamped  Isotropic  Arch 


As  can  be  seen,  the  initial  load-up  for  each 
formulation  is  approximately  the  same.  The  SHELL  code 
calculates  a  higher  buckling  load  because  it  includes 
nonlinear  inplane  displacement  terms  in  the  strain  relations 
that  are  not  included  by  Belytschko  and  Glaum  which  also 
does  not  include  through  the  thickness  shear  terms.  As  the 
displacement  becomes  large  a  greater  effect  of  extensibility 
in  the  SHELL  element  becomes  present. 
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Dynamic  Response  of  Isotropic  Cylindrical  Shell 


A  three  inch  thick  isotropic  cylindrical  shell  (shown 
in  Fig  3-7)  was  subjected  to  a  uniformly  distributed  half 
sine  wave  impulsive  loading  with  peak  intensity  of  90  psf  as 
shown  in  Fig  (3-8) . 
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Figure  3-8.  Loading  for 

Isotropic  Shell 


Results  are  compared  to  existing  results  by  Clough  and 
Wilson  (5) ,  where  through  the  thickness  parabolic  shear 
strain  is  not  included.  The  two  straight  longitudinal  edges 
were  assumed  free  and  the  two  circular  edges  were  assumed  to 
be  supported  on  diaphragms.  Material  properties  of  E  - 
10.0E7  and  v  =  0  are  used.  As  done  previously,  1/4  panel 
symmetry  is  taken  advantage  of.  The  time  increment  in  each 
load  step  was  0.025  sec  which  is  approximately  1/25  the 
first  natural  period  of  the  shell. 

The  initial  two  analyses  were  accomplished  without 
considering  the  inertia  force,  i.e.  setting  the  mass  terms 
to  zero,  thus  they  were  effectively  "static”  analyses  even 
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though  the  displacement  changed  with  time.  Both  geometric 
linearity  and  non-linearity  were  considered  by  using  SHELL’S 
capability  to  either  include  or  exclude  the  nonlinear  terms 
in  the  strain-displacement  relations.  The  results  show  good 
agreement  and  are  presented  in  Fig  (3-9) . 


Figure  3-9.  Displacement  of  Pt  A  vs  Time,  ’’Static” 


3-10 


I 

|  Two  additional  analyses  were  carried  out  in  which  the  mass 

density  of  the  material  is  considered.  The  dynamic 
*  deflections  are  larger  than  static,  as  would  be  expected. 

|  Comparison  with  results  from  Clough  and  Wilson  are  shown  in 

Fig  (3-10) .  The  agreement  of  results  is  due  to  the  fact 
j  that  transverse  shear  deformations  are  small  for  a  shallow 

thin  isotropic  arch  and  the  loading  is  well  below  the 
I  critical  load  which  would  cause  snap-through. 


|  Figure  3-10.  Displacement  of  Pt  A  vs  Time,  "Dynamic" 

I 

I 
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Critical  Time  Step 


Before  any  further  results  are  presented,  it  is 
important  to  discuss  a  parameter  used  in  the  numerical 
analysis  of  non-linear  dynamics  called  the  critical  time 
step  (At) .  Shell  structures  usually  result  in  a  so-called 
"stiff"  system  of  differential  equations  where  the  highest 
eigenvalues  determine  the  numerical  stability  of  the 
explicit  integration  method  (1).  Therefore,  a  very  small 
time  step  must  be  used.  If  the  chosen  At  exceeds  the 
critical  value,  a  small  error  will  be  magnified  with  each 
time  step  and  computed  displacement  will  grow  very  rapidly 
and  a  numerical  problem  will  be  noticed  after  only  a  few 
time  steps.  If  At  is  less  than  the  critical  value  an  error 
will  not  accumulate  and  the  solution  obtained  will  be 
reasonably  accurate  to  the  degree  of  tolerance  used  in  the 
numerical  scheme.  Katona  and  Zienkiewkz  (12)  suggest  using 
a  At  at  most  1/12  times  the  first  natural  frequency  of  the 
structure.  Reference  (1)  proposes  several  different  methods 
for  estimating  the  critical  time  step  but  concludes  that 
once  an  initial  estimate  is  made  trial  and  error  will  most 
likely  be  necessary  to  further  refine  the  critical  time 
step.  Dinkier  and  Kroplin  (9)  also  point  out  that  errors 
will  grow  uncontrollable  unless  a  small  time  step  is  used 
but  they  also  do  not  offer  any  explicit  methods  for 
determining  At. 
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Laminated  Arch 


In  this  section,  both  the  static  and  dynamic  analysis 
of  a  simple  supported  laminated  arch  subjected  to  a  center 
point  load  are  presented. 

Material .  The  analysis  is  performed  on  laminate 
constructed  of  high  strength/low  modulus  Hercules  AS4-3501-6 
graphite  epoxy  with  the  following  ply  physical  properties 
and  dimensions: 

ply  thickness  =  0.005  in. 

E  1  =  18.844E6  psi 

E  2  =  E  3  =  1.468E6  psi 

G  i2  =  G  13  =  0.91E6  psi 

G  23  =  0.45E6  psi 

v  12  =  0.28 

V21  =  V31  =  °*0218 

mass  density  =  0.00015088  slugs/ in  3 

Arch  Geometry.  A  ply  lay  up  of  [0  ^90  6]  s  is  used 
which  has  24  plies  and  results  in  a  total  thickness  of  .12 
inch.  As  shown  in  Fig  (3-11)  a  radius  of  curvature  of  12 
in.  along  with  an  open  angle  0  of  1  radian  is  used.  Width 
is  the  same  order  of  thickness  at  .12  in.  Mesh  size  for  the 
quarter-arch  ranged  from  1”  x  .06”  to  .25”  x  .06”  as  shown 
in  Fig  (3-12) . 
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Figure  3-11.  Laminated  Arch 


Static  Analysis.  The  SHELL  code  was  used  to 


investigate  the  static  snap-through  of  the  laminated  arch. 
The  static  load-deflection  curve  is  illustrated  in  Fig  (3- 
13) ,  which  shows  that  the  critical  snapping  load  is 
approximately  11.5  lb. 


Figure  3-13.  Load  vs  Displacement  of  Center  Node, 
Laminated  Arch 


If  a  load  slightly  greater  than  11.5  lb  were  applied  it 
can  be  concluded  from  the  curve  that  the  shell  should 
collapse  and  the  center  node  should  displace  approximately 
2.7  in.  Figure  3-14  shows  a  profile  of  the  centerline  of 
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the  arch  at  pt  A  and  pt  B  from  Fig  (3-13)  along  with  the 
original  configuration. 


The  arch  instantaneously  snaps  through  from  point  A  to  B  and 
is  unstable  until  reaching  point  B  where  it  is  able  to 
resume  sustaining  a  load.  As  can  be  seen  from  Fig  (3-14) 
the  arch  at  point  A  is  partly  above  and  below  its  original 
configuration.  This  suggests  that  the  part  above  is  in 
compression  and  the  part  below  is  in  tension.  Once  the  arch 
has  fully  collapsed  it  is  entirely  in  tension.  Each  point 
is  below  the  original  configuration. 
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Dynamic  Analysis.  The  same  problem  subjected  to 
dynamic  loading  is  investigated.  In  the  first  case 
investigated,  the  applied  dynamic  loading  history  is  shown 
in  Fig  (3-15) ,  where  the  maximum  load  is  smaller  than  the 
static  critical  snapping  load. 


Figure  3-15.  Loading  History 
for  Laminated  Arch,  Pre¬ 
buckling  Load 


The  response  is  a  regular  smooth  periodic  vibration  as  shown 
in  Fig  (3-16) .  One  curve  shows  a  case  where  no  damping  is 
allowed  and  thus  the  amplitude  of  each  period  is  the  same 
and  the  other  case  shows  that  if  damping  is  artificially 
imposed  the  steady  state  solution  will  result.  A  time 
increment  of  .0002  sec  was  used.  Since  the  displacement 
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versus  time  curve  displays  good  periodicity  and  the  loading 
is  less  than  critical  where  non-linear  effects  are  not 
dominant,  this  time  step  was  judged  to  be  acceptable.  As 
can  be  seen  from  the  undamped  case,  approximately  17  points 
are  plotted  for  each  period  which  satisfies  Katona  and 
Zienkiewicz'  (12)  criteria  of  at  least  12  per  period.  Thus, 


TIME  (sec) 


Figure  3-16.  Displacement  vs  Time  of  Center  Node 
for  Laminated  Arch,  Pre-buckling  Load 


the  natural  frequency  can  be  approximated  at  294  Hz.  An 
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eigenvalue  analysis  done  by  the  SHELL  code  calculates  the 
natural  frequency  at  499  Hz.  Since  the  eigenvalue  analysis 
uses  only  the  linear  terms  in  the  stiffness  matrix  one  would 
expect  any  calculations  using  the  non-linear  terms  to 
predict  a  lower  natural  frequency. 

In  the  second  case  analyzed,  a  maximum  load  of  12  lb  is 
applied  as  shown  in  Fig  (3-17) ,  which  is  slightly  greater 
than  the  critical  snapping  load  of  11.5  lb. 


Figure  3-17.  Loading  History 
for  Laminated  Arch,  Post- 
buckling  Load 


A  time  increment  of  0.0002  sec  is  again  used.  During 
this  phase  of  the  research  it  was  discovered  that  the 
percent  convergence  tolerance  of  the  numerical  technique  is 
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very  critical  when  investigating  dynamic  behavior  due  to 
loads  above  critical.  Tolerances  of  0.1  -  0.2  percent  are 
acceptable  when  the  load  is  sub-critical  but  this  tolerance 
range  causes  errors  to  grow  uncontrollable  when  the  load  is 
above  critical.  It  becomes  obvious  that  the  solution 
"breaks  down"  before  representing  a  realistic  dynamic 
pattern.  Fig  (3-18)  shows  a  comparison  of  the  time  versus 
center  node  displacement  history  obtained  using  a  percent 
convergence  tolerance  of  0.1  and  0.0001.  The  0.1  percent 
curve  "breaks  down"  shortly  after  reaching  the  peak 
displacement.  As  can  be  seen  error  grows  uncontrollably  and 
the  curve  no  longer  represents  a  realistic  solution.  The 
0.0001  percent  curve  calculates  two  realistic  vibration 
periods  before  it  was  no  longer  able  to  converge  to  the 
desired  tolerance. 
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I 

I 


TIME  (sec) 


Figure  3-18.  Center  Node  Displacement  vs  Time, 
Laminated  Arch 


Fig  (3-19)  shows  the  time  vs  displacement  history  for 
both  an  undamped  and  artificially  damped  case  using  the 
percent  convergence  tolerance  of  0.0001.  The  damping 
coefficient  imposed  is  0.1.  Because  the  non-linear  parts  of 
the  strain-displacement  relations  are  dominant  and  the 
solution  technique  is  numerical,  a  pure  sinusoidal  curve  is 
not  seen  but  it  is  approximated.  The  reason  as  to  why  the 
two  cases  do  not  average  the  same  displacement  is  as  yet 
unresolved.  As  can  also  be  seen,  the  initial  curve  leading 
up  to  the  maximum  displacement  is  concave  up.  This 
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demonstrates  that  the  nonlinear  kinematics  are  producing  a 
more  flexible  arch  response  relative  to  linear  relations. 


TIME  (sec) 


Figure  3-19.  Displacement  vs  Time  of  Center  Node, 
Laminated  Arch,  Post-buckling  Load 


Fig  (3-20)  shows  a  profile  of  the  arch  centerline  at  points 
A,  B,  and  C  from  Fig  (3-19) .  As  can  be  seen,  the  centerline 
at  point  B,  which  is  at  approximately  the  average  amplitude, 
closely  resembles  that  from  Fig  (3-14) .  However,  it  is 
unclear  why  the  displacement  at  point  B,  which  is  at  the 
average  amplitude  for  the  undamped  case,  is  2.5  inches  while 
the  static  curve  suggests  a  collapse  displacement  of  2.7 
inches. 
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Figure  3-20.  Profile  of  Laminated  Arch  Centerline 

Fig  (3-20)  also  shows  that  only  the  middle  third  of  the  arch 
appears  to  be  in  vibration.  There  are  two  counter  f lecture 
points  approximately  1/3  in  from  each  boundary  where  the 
bending  moment  remains  near  zero.  Also  shown  is  a  profile 
of  the  arch  corresponding  to  point  X  on  Fig  (3-19)  which  is 
immediately  after  the  onset  of  the  maximum  load.  This 
profile  closely  resembles  the  pre-buckling  profile  from  Fig 
(3-14) .  This  profile  does  not  appear  immediately  after  the 
onset  of  the  critical  load  because  the  dynamic  calculations 
take  into  account  inertia  forces  but  it  does  indicate  the 
start  of  the  snapping  phenomena. 


Laminated  Cylindrical  Shell 

The  final  case  investigated  is  that  of  a  laminated 
cylindrical  shell  subject  to  a  center  point  load.  The 
material  studied  is  the  same  as  the  laminated  arch  and 
presented  again  below: 

ply  thickness  =  0.005  in. 

E1  =  18.844  x  106  psi 
E2  =  E3  =  1.468  x  106  psi 

G12  =  G13  =  0>91  x  1q6  Psi 
G23  =  0.45  x  106  psi 
v12  =  0.26 

V21  =  V31  =  °-0218 

mass  density  =  0.00015088  slug/in3 

Geometry  and  Boundary  Conditions.  The  same  ply  lay-up 
as  the  arch  was  used  of  [06/906]s  resulting  in  an  overall 
shell  thickness  of  0.12  in.  Again  a  radius  of  curvature  of 
12  inches  is  used  along  with  an  open  angle  0  of  1  radian  as 
shown  in  Fig  (3-21) .  A  length  of  18  inches  is  used  so 
results  can  be  compared  to  Silva  (23)  who  ran  the  same 
static  analysis  on  the  SHELL  code.  Both  straight  edges  were 
allowed  to  be  simply  supported  while  both  curved  edges  were 
allowed  to  be  free.  This  combination  of  geometry,  ply- 
layup,  and  boundary  conditions  were  chosen  because  it  was 
shown  by  Silva  to  fully  collapse.  "Full"  collapse  is  said 
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to  occur  when  the  free  ends  also  turn  under  which  implies 
that  every  point  in  the  shell  has  displaced  below  its 
original  configuration. 


Finite  Element  Mesh.  As  seen  in  Fig  (3-21) ,  Silva  uses  an 
88  element  1/4  panel  mesh  to  achieve  a  static  load  versus 
displacement  curve.  In  an  effort  to  reduce  the  necessary 
computer  run  time,  particularly  for  the  dynamic  analysis, 
the  present  work  reduces  the  number  of  elements  to  a  49 
element  mesh.  Although  the  number  of  degrees  of  freedom  in 
this  mesh  does  not  fall  within  the  convergence  criteria  set 
by  Silva  (23) ,  it  does  fall  within  5  percent  of  his  criteria 
based  on  peak  load.  Silva  uses  elements  ranging  in  size 


from  l"xl"  to  V'xV  while  the  present  work  uses  elements 
that  range  in  size  from  2"xl"  to  V'x^" .  It  should  be  noted 
that  in  the  area  of  point  load  application  the  element  size 
for  both  meshes  is  V'xV*. 
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Figure  3-22.  Mesh  Comparison  for  Laminated  Shell 
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Static  Analysis.  Fig  (3-23)  shows  the  static  load 
versus  displacement  curves  for  both  Silva's  analysis  and  the 
present  work.  The  initial  load  up  and  critical  load  are 
approximately  the  same  but  with  less  finite  elements  the 
post  buckling  behavior  is  slightly  more  difficult  to  track. 
The  critical  snapping  load  is  calculated  to  be  approximately 
2800  lb.  The  initial  post  buckling  deflection  of  the  center 


Figure  3-23.  Load  vs  Displacement,  Cylindrical  Shell 
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Dynamic  Analysis.  The  same  problem  subjected  to 
dynamic  loading  is  investigated.  For  the  first  analysis, 
the  applied  dynamic  loading  history,  shown  in  Fig  (3-23) ,  is 
smaller  than  the  static  critical  snapping  load.  The 
response,  as  shown  in  Fig  (3-24),  is  periodic,  but  not 


Figure  3-24.  Loading  History 
for  Cylindrical  Shell, 
Pre-critical  Load 


purely  sinusoidal.  The  fact  that  the  curve  is  not  purely 
sinusoidal  can  perhaps  be  accounted  for  due  to  the  large 
number  of  calculations  performed,  the  non-linearity  of  the 
strain-displacement  equations,  the  precision  of  the  computer 
code  (double  precision  is  used)  and  the  round-off  error  of 
the  computer.  The  average  displacement  for  each  period  does 
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however,  match  that  for  the  static  load  vs  deflection  curve 
of  Fig  (3-23)  which  indicates  that  steady  state  has  been 
reached. 

A  time  step  of  0.0002  seconds  is  presented  however  the 
results  were  the  same  when  a  time  step  of  0.0001  sec  was 
used  thus  convergence  is  assumed.  By  approximating  the  time 
between  peaks,  a  natural  frequency  of  225  Hz  is  estimated. 
Again  an  eigen  value  analysis  predicts  a  higher  natural 
frequency  of  505  Hz. 


Figure  3-25.  Displacement  vs  Time  of  Center  Node, 
Cylindrical  Shell,  Pre-critical  Loading 
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In  the  second  analysis,  the  applied  dynamic  loadinj  history 
is  greater  than  the  static  critical  snapping  load  as  shown 
in  Fig  (3-26).  Time  steps  of  0.0002  sec  and  0.0001  sec  were 
used.  A  percent  tolerance  of  0.001  was  imposed. 


Figure  3-26.  Load  History  for 
Cylindrical  Shell,  Post- 
Buckling  Load 


The  vibration  response  of  the  shell  is  shown  in  Fig  (3-27) . 
The  curve  representing  the  time  step  of  0.0002  sec  "breaks 
down"  after  reaching  the  maximum  displacement  and  is  not 
able  to  show  any  reasonable  periodicity.  By  cutting  the 
time  step  in  half  to  0.0001  seconds  the  curve  is  able  to 
show  approximately  one  and  one-half  cycles  before  it  also 
reaches  a  point  where  the  error  goes  uncontrollable. 
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However,  during  the  one  and  one-half  periods  the  average 
amplitude  of  2.7  inches  matches  that  of  the  static  load 
displacement  curve  of  Fig  (3-23) . 


Figure  3-27.  Displacement  vs  Time  of  Center  Node, 
Cylindrical  Shell,  Post-buckling  Load 


At  this  point  it  is  as  yet  unresolved  as  to  whether 
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decreasing  the  time  step  even  further  would  yield  even  more 
realistic  results  or  the  limits  of  the  numerical  precision 
of  the  computer  have  been  reached. 
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IV.  CONCLUSIONS /RECOMMENDATIONS 


Conclusions 

The  instability  of  thin  laminated  cylindrical  shells 
acted  upon  by  transverse  loads  is  clearly  a  dynamic 
phenomena  for  which  analytical  prediction  and  physical 
reality  have  yet  to  be  fully  reconciled.  Some  authors 
prefer  to  investigate  the  post  buckling  state  using  a  static 
analysis  while  others  argue  that  only  a  dynamic  analysis  is 
able  to  cope  with  the  difficult  phenomena  of  transient 
multimode  buckling  between  the  pre-  and  post-buckling  range. 
It  is  apparent  that  this  argument  becomes  very  significant 
when  instability  of  a  shell  is  considered  with  transverse 
loading.  It  is  also  a  major  characteristic  when  one 
investigates  cylindrical  panels  acting  under  axial 
compression.  With  the  SHELL  code,  it  is  possible  to  use  the 
same  set  of  constitutive  laws,  strain-displacement 
relations,  and  finite  element  formulations  to  do  both  a 
static  and  dynamic  analysis  on  a  particular  shell.  The 
resulting  comparison  gives  insight  into  the  post  buckling 
behavior  of  shells  and  their  ability  to  sustain  a  load  after 
collapse.  An  exclusively  static  analysis  is  a  useful  tool 
for  predicting  the  immediate  post  buckling  state  in  a  dead 
load  situation.  However,  such  an  analysis  considers  only 
static  equilibrium  states.  Beyond  the  critical  load  point 
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it  is  concerned  with  equilibrium  states  which  may  or  may  not 
occur  as  the  shell  dynamically  seeks  another  configuration 
in  which  it  can  resume  supporting  a  load.  On  the  other 
hand,  a  dynamic  analysis  is  not  concerned  with  these 
artificial  equilibrium  states  and  can  trace  the  behavior  of 
the  shell  all  the  way  through  the  buckling  process.  It 
should  be  noted  though,  that  the  dynamic  analysis  requires 
significantly  more  computer  time  than  a  static  analysis.  In 
a  numerical  scheme,  the  time  step  must  be  small  enough  to 
avoid  numerical  instability.  It  has  been  shown  that  in 
using  the  SHELL  code,  the  response  of  a  shell  resulting  from 
a  dynamic  steady  state  analysis  subjected  to  a  step  load, 
matches  the  displacement  on  the  snap  through  load  versus 
displacement  curve. 

Recommendations 

There  are  many  parameters  involving  the  use  of  the 
SHELL  dynamic  code  which  lend  themselves  as  areas  of  further 
study.  The  effect  of  boundary  conditions,  other  than  simply 
supported,  require  further  study.  For  example,  a  shell  that 
is  clamped  on  the  edges  will  respond  differently  and  is 
perhaps  more  common  in  aircraft  structures.  Also,  various 
load  conditions  such  as  line  and  distributed  loads  should  be 
investigated.  An  impulse  load  can  also  be  used  other  then  a 
step  load.  The  use  of  a  damping  factor  requires  further 
study  particularly  in  trying  to  get  a  physical  meaning  of 
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the  phenomena.  Perhaps  the  most  important  area  requiring 
further  study  is  that  of  the  numerical  integration  scheme 
and  its  inherent  numerical  instability.  More  investigation 
into  the  tolerance  and  the  critical  time  step  are  needed. 
Also,  the  use  of  different  parameters  for  the  beta-m 
integration  technique  might  cause  convergence  to  a  solution 
at  a  faster  rate. 
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